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Colloidal particles that are confined to an interface such as the air-water interface are an example 
of a two-dimensional fluid. Such dispersions have been observed to spontaneously form cluster and 
stripe morphologies in certain systems with isotropic pair potentials between the particles, due to 
the fact that the pair interaction between the colloids has competing attraction and repulsion over 
different length scales. Here we present a simple density functional theory for a model of such a two- 
dimensional fluid. The theory predicts a bulk phase diagram exhibiting cluster, stripe and bubble 
modulated phases, in addition to homogeneous fluid phases. Comparing with simulation results 
for this model from the literature, we find that the theory is qualitatively reliable. The model 
allows for a detailed investigation of the structure of the fluid and we are able to obtain simple 
approximate expressions for the static structure factor and for the length scale characterising the 
modulations in the microphase separated phases. We also investigate the behaviour of the system 
under confinement between two parallel hard walls. We find that the confined fluid phase behaviour 
can be rather complex. 
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I. INTRODUCTION 

Many soft matter systems are capable of displaying 
spatially modulated structures. Well known systems such 
as solutions of amphiphilcs and diblock copolymers can 
exhibit several microphase ordered phases p~j. This self 
organisation originates from competing interactions oc- 
curring over different length scales and particle orienta- 
tions. Systems of colloidal particles confined to the air- 
water interface are also capable of exhibiting microphase 
separation [HO [3]. In these two dimensional fluids, the 
self organisation is driven by a competition between short 
ranged attractive interactions, which gather the particles 
together, but at longer ranges there is a repulsive inter- 
action between the particles which breaks the system up 
into high density islands of particles with low density re- 
gions in between the islands. The shape of these islands 
depends on the surface coverage (density) . At lower den- 
sities the particles form an array of (roughly) circular 
clusters, but at higher densities, they form a series of 
elongated parallel stripes [3j [4J . This self organisation is 
rather striking, when bearing in mind that the interac- 
tions between the individual pairs of particles is only a 
function of the particle separation and not of any orien- 
tational degrees of freedom, which play an important role 
in the microphase ordering of, for example, amphiphilic 
systems. 

Suspensions of colloidal particles (i.e. three dimen- 
sional fluids) with effective pair interactions between the 
particles of a similar form have also been synthesised 
[5j[6jl3[8]. These systems also display various microphase 
separated phases including cluster phases. An overview 
of the various studies of models for these three dimen- 
sional systems is outlined in Ref. [§]• 

The focus of this paper is a model two-dimensional 



fluid, first proposed in Ref. [3], and studied in detail via 
Monte Carlo computer simulations in a number of pub- 
lications by Imperio and Reatto pH HU [El [13]. This 
system is composed of particles with a hard core of di- 
ameter er, interacting via the following pair potential: 

v(r) = < ( . (1) 
I w(r) r > a, 

where 

w(r) = ~^ exp (~£) + ^ exp (-i)- (2) 

The parameter subscripts a and r stand for "attraction" 
and "repulsion" respectively. In much of the work pre- 
sented here we will follow Imperio and Reatto and se- 
lect the following values for the pair potential parame- 
ters: R a = Iff, R r = 2cr, and e a = e r = e. Thus e 
will be the unit of energy in the system. These parame- 
ters are chosen so that the integral over w(r) is zero, i.e. 
a = — J dr w(r) = 0. 

The studies in Refs. [U [lOj HU [12l H3] found a rich 
phase behaviour displayed by this model fluid, which mir- 
rored the behaviour observed in the experiments [31 UJ ■ 
At high temperatures, the contribution from w(r) to v{r) 
is negligible and so the fluid properties are wholly deter- 
mined by the hard-core part of the potential v(r) and 
the system is effectively just a fluid of hard-disks. On 
lowering the temperature, however, one finds that on in- 
creasing the fluid density at fixed temperature, there is 
a sequence of phase transitions. At low densities, the 
particles are uniformly distributed across the surface; we 
denote this the vapour phase. As the density is increased 
we find a transition to the cluster phase, in which the 
particles gather to form circular islands. The number of 
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particles in each island may vary from just a few to as 
many as a hundred or more, depending on the particular 
values of the pair potential parameters in Eq. The 
clusters are equilibrium structures and the probability of 
a particle moving from one cluster to another is > 0. 
There is no long range order between the clusters - the 
system forms a fluid of clusters. However, as the den- 
sity is further increased, the clusters 'freeze' so that the 
clusters arrange themselves onto a triangular superlattice 
[TT] . On further increasing the density of particles on the 
surface, there is a transition to a phase in which the par- 
ticles form parallel stripes [H [TUJ E] • Further increasing 
the density, finds the system exhibiting a 'bubble' phase, 
in which there are low density voids amongst the parti- 
cles on the surface. At even higher densities the system 
forms a uniform liquid phase |14j . It is worth emphasis- 
ing again that for some temperatures, these non-uniform 
structures are equilibrium (ergodic) states, having a non- 
zero probability for a particle to move from one cluster 
or stripe to another. However, as the temperature is 
decreased, this transition probability decreases. At low 
temperatures, one finds the same sequence of phase tran- 
sitions, but for these temperatures the particles can be 
frozen within the clusters or stripes. The possibility of 
a related glass transition should not be ruled out either 

It should also be recalled that pattern formation in 
two dimensional model fluids composed of particles in- 
teracting via isotropic pair potentials is not restricted to 
generalisations of Eqs. (JlJ and ([2| - i.e. potentials with 
a short ranged attraction and a longer ranged repulsion. 
Striped patterns have also been observed in systems of 
particles interacting via a pair potential having a hard- 
core plus a longer ranged purely repulsive 'shoulder' po- 
tential [15 and also in these systems when an additional 
attractive contribution is added to the potential, beyond 
the repulsive shoulder - i.e. potentials with a short ranged 
repulsion and a longer ranged attraction [16] . 

In this paper we develop and apply a simple density 
functional theory (DFT) for the model fluid defined by 
Eqs. <Q and @. In Sec.|H] we describe the DFT. In Sec. 
|III| we use the DFT we investigate the structure of the 
fluid and, in particular, we obtain a simple expression 
for the static structure factor and also for the character- 
istic length scale associated with the modulations in the 
non-uniform phases. In Sec. |IV| we calculate the phase 
diagram predicted by the simple DFT and find that it 
exhibits a cluster, stripe and bubble phase, which is in 
qualitative agreement with the simulation results of Refs. 
[IU1[II1|T|. In Sec. [V] we use the DFT to study the 
phase behaviour of the fluid when confined between two 
parallel hard walls. Depending on the fluid density and 
temperature, we find that as the separation between the 
two walls is varied, the equilibrium density profile can 
vary significantly due to the need for the length scale of 
the density modulations in the fluid to be commensurate 
with the distance between the walls. Finally, in Sec. |VI| 
we discuss our results and draw some conclusions. 



II. A DFT FOR THE SYSTEM 

The structural and thermodynamic properties of the 
system may be obtained from the grand potential func- 
tional Gil HH1: 



n[p(v)} = T[p{v)] + J drp(r)(F ext (r) - p), (3) 

where p(r) is the fluid one body density, V ext (r) is the ex- 
ternal potential, p is the chemical potential and .F[/>(r)] is 
the intrinsic Hclmholtz free energy functional. The grand 
potential of the system for a given p and V ext (r) , denoted 
f2, is the minimal value of f2[/>(r)], and the equilibrium 
one-body density profile p(r) is that which minimises the 
grand potential functional QTJ [18]. As usual, ,F[/o(r)] 
is an unknown functional and we assume the following 
simple mean-field approximation for this quantity: 

T[p(r)] = J drp(r)/(p(r)) 

+ l -JdvJdv'p{v)p{v')w{v-v l ) (4) 

where f(p) is the Helmholtz free energy per particle of a 
uniform fluid of hard-disks with bulk density p. This local 
density approximation (LDA) for the reference hard disk 
contribution includes the (exact) ideal gas contribution 



^d[p(r)] = k B T / drp(r)(ln[A 2 p(r)] - 1), 



(5) 



where A is the thermal De Brogue wavelength. For sim- 
plicity, we will use the scaled particle approximation, 

f3f(p) = ln(A 2 p) - 2 - ln(l - ij) + (1 - (6) 

where /3 = 1/fcgT is the inverse temperature and 77 = 
7t/5ct 2 /4 is the packing fraction. There exist better ap- 
proximations than Eq. (|6|, such as that obtained in Ref. 
[19] , but for the densities of interest here, the more sim- 
ple expression in Eq. (|6| is sufficiently reliable. The mean 
field contribution to the free energy in Eq. Q, from the 
tail of the pair potential, is justified on the basis that 
w(r) is fairly long ranged and slowly varying. Note that 
Eq. Q does not provide a reliable approximation for the 
Helmholtz free energy in cases when the fluid density p(r) 
varies strongly on length scales ~ a, i.e. for describing ef- 
fects from packing of the hard cores of the particles. To 
take account of such correlations, one must implement a 
non-local approximation for the reference hard-disk free 
energy functional, along the lines of that used in Ref. [2D] 
for the three-dimensional counterpart of the present sys- 
tem. However, in the present study, where our interest is 
in the cluster, stripe or bubble structures, for which the 
fluid density profile p(r) varies over length scales 3> <r, 
Eq. Q suffices [2T]. 
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III. STRUCTURE OF THE UNIFORM FLUID 



The direct pair correlation function, denned as follows 



c(r,r') = -/3 



5\F[p{v)]-T ld [p{*)]) 
Sp(r)6p(r>) 



(7) 



can be used to characterise two-point correlations in the 
inhomogeneous fluid. In the homogeneous bulk fluid, 
where p(r) = p, we find c(r,r') = c(|r — r'|) = c(r). More 
conventionally one characterises the fluid structure by ei- 
ther considering the radial distribution function g(r), ob- 
tained from c(r) via the Ornstein Zernike equation [17], 
or the static structure factor 



S(k) 



1 - pc(k) ' 



(8) 



where c(k) is the Fourier transform of c(r). 

For the present model, Eqs. Q - Q together yield 
the following expression for the pair direct correlation 
function 



c(r.r') = -/3 



2/'(p(r)) + p(r)/"(p(r))- 



k B T 
P(r) 



-/3u>(r — r'), 



5{v-v') 
(9) 



where S(r) is the Dirac delta function and /' and /" 
are the first and second derivatives of / with respect to 
p. Combining Eqs. (|6| - ^ we obtain the following 
expression for the static structure factor: 



S(k) 



1 



a(rj) + p(3w{k) ' 



where a{ij) = (1 + 77) /(l — ij) 3 and where 



w(k) 



2ne a a 2 



2ire r a 2 



(1 + Rlk 2 ) 3 / 2 (l + i?2fc 2 )3/2- 



(10) 



(11) 



Due to the local density approximation for the hard-disk 
contribution to the free energy functional ^ , this expres- 
sion for S(k) is strictly only applicable for wave numbers 
k < 2ir/a. In Ref. [TO], the authors showed that the 
RPA approximation for S{k) is fairly reliable for describ- 
ing the fluid structure. The RPA consists of making the 
approximation c(fc) = Chd(k) — (3w{k), where Chd{k) is 
the Fourier transform of the hard disk pair direct corre- 
lation function, c^lV). For the present model fluid, the 
RPA approximation yields a structure factor which is re- 
liable for all wave numbers k. The result in Eq. ( 10 ) only 
coincides with the RPA result for small wave numbers 
k < 2ixja. Thus the present theory effectively only in- 
cludes the inter-cluster correlation contributions to S(k) 
and does not include the intra-cluster correlation contri- 
butions. 

As an aside, we now briefly consider the structure and 
phase behaviour of a system for which e r = 0, i.e. when 
the long range repulsive contribution is no longer present 



and (-a) = J drw(r) = w{k = 0) = -2-Kt a a 2 < 0. 
In this case, the system exhibits only bulk vapour- liquid 
phase coexistence and no microphase separation. Owing 
to our simple approximation for the free energy Q, it 
is straightforward to determine the phase diagram and 
calculate the binodal and spinodal. For example, the 
spinodal corresponds to the divergence S(k — 0) — > 00; 
i.e. this is when the denominator in Eq. ( |10| equals zero 
for k — 0, giving the following expression for the spinodal 
temperature as a function of density: 



k B T _ p{l-r,f 



(1 + ??) 



(12) 



The critical point corresponds to the maximum of this 
curve, i.e. when — ?7) 3 /(l + if)] = 0. This yields a 

critical packing fraction r\ c = (V7 — 2)/3 ~ 0.274. Note 
that one could also obtain Eq. ( fl2"| starting directly from 
the free energy functional Q. In bulk, this gives the 
following expression for the bulk free energy per particle: 
ftot(p) = f(p) ~ pot/2. From this, we may obtain the 
pressure in the system: 



P 



P 



fdft, 



V dp 
pk B T 



p 2 a 



(13) 



The spinodal corresponds to the locus (dP/dp)x = 0, so 
taking a further derivative we obtain Eq. ( |12[ ). Such con- 
sistency between the structural and the thermodynamic 
routes to the phase behaviour comes as a consequence 
of deriving all quantities from a free energy functional 
and is one of the advantages of using the DFT approach 
rather than integral-equation theory based approaches 

mi [m cia. 

We return now to the more general case when e r =/= 0. 
First, we focus on the small wave number k behaviour of 
S(k). Taylor expanding S(k) in Eq. (10 1 around k = 
we find: 



S(k) = S - jSf t k 2 + 0{k A ) 



(14) 



where So = S(k = 0) = [a{rf) — p^a]^ 1 and 7 = 
37TCT 2 pf3(e a R 2 - e r Rl). When 7 > 0, i.e. when (e a R 2 - 
e r R 2 ,) > 0, then S(k = 0) is a maximum and the phase 
behaviour of the system is the same as in the case when 
e r = 0, provided a > 0, where a = — J drw(r) — 
2-Ko 2 (t a — e r ), and the spinodal curve is given by Eq. 
(12 1. However, when 7 < 0, i.e. when (e a R 2 a — e r R 2 ) < 0, 
then S(k = 0) is a minimum and there is a maximum in 
S(k) at k — k c ^ 0. This peak in S(k) at k c indicates 
a propensity in the fluid for having density modulations 
with wavelength 2ir/k c . In fact, in certain portions of 
the phase diagram, the fluid may become unstable with 
respect to density fluctuations of wavelength 2n/k c and 
it then follows that the system exhibits microphase sep- 
aration to the various different modulated fluid phases, 
in this portion of the phase diagram. This instability is 



4 



indicated by the divergence of the peak in S(k) at k c . 
We denote the locus in the phase diagram at which this 
occurs the A-line - i.e. the A-line is defined as the locus 
where S(k = k c ) oo [HI [HJ [231 [Ml HSJ [26]. 

Within the present DFT we may obtain an expression 
for k c and also for the temperature on the A-line, as a 
function of the density, as follows: the wave number k c 
is defined as the value of k where S(k) is a maximum, 
i.e. where the quantity [a{rf) + pj3w{k)\ is a minimum 
[recall Eq. ( 10 1] . This occurs when -§r[a(i]) + pj3w{k)] = 



0.28 



pP-^w(k) = 0. Solving this equation we find that there 
are two solutions. The first at k = 0, corresponding to 
the minimum in S(k) and the second at k = k c , where 



k r . — 



' r-i 
m - rm ' 



(15) 



and where T = (e r R^/e a Ra) 2 ^ 5 ■ For the set of pair poten- 
tial parameters that we focus on here, R a = la, R r = 2a, 
and e a = e r = e, we find k c a ~ 0.573 so that the length- 
scale characterising the density modulations in the mi- 
crophases = 2ir/k c ~ Her, as pointed out in Ref. [ID] . 



Comparing the result in Eq. (15 1, with the results from 



Monte-Carlo computer simulations, for various different 
sets of choices for the pair potential parameters, one finds 
that Eq. (151 is generally reliable and in good agreement 



with the simulation results 

The A-line corresponds to the locus where [a(rf) + 
pj3w{k — k c )\ — 0. Rearranging this we obtain the fol- 
lowing expression for the temperature along the A-line: 



k B T pil-v? 



(1 + V) 



(16) 



where A = —w(k = k c ) > 0. Note the strong similarity 
between this expression and Eq. ( 12 ) for the spinodal. 
For the pair potential parameters R a = ler, R r = 2a, 
and e a = e r = e, the A-line is displayed in Fig.[l] See also 
Fig. 14 in Ref. [TDJ, where the A-line from the full RPA 
is displayed. Inside the A-line, the homogeneous fluid is 
unstable with respect to density fluctuations with wave 
number k c . The nature of the phases that occur in this 
region of the phase diagram is the topic of the following 
section. 



IV. BULK PHASE BEHAVIOUR 

For a given state point, the equilibrium fluid density 
profile is that which minimises the grand potential (|3]); 
i.e. it is the solution to the following equation [171 H8jT 



6Sl[p(r)} 



= 0. 



(17) 



f(p(T)) + P (r)f(p(r))+ I dr'p(r')w(T-r')+V ext (v) = p. 

(18) 



Sp(r) 

From Eqs. Q and ([TFl we obtain 




FIG. 1: (Color online) Phase diagram in the density- 
temperature plane. The shaded regions denote two-phase co- 
existence regions between the various phases and the dashed 
line is the A-line. 



We discretise the density profile p(r) on a two- 
dimensional Cartesian grid and, for a given value of the 
chemical potential p, we may then solve Eq. (18) for the 
fluid density profile, by using a simple iterative numer- 
ical procedure. For the bulk case, where the external 
potential V ex t(r) — 0, we solve using periodic boundary 
conditions in both Cartesian directions. Since in this 
case the uniform density profile p(r) = p is always a 
solution of Eq. ( 18 1 , to find the density profile for the 
modulated phases we investigated three approaches: (i) 



choose an initial guess for the input density profile with a 
cluster/stripe/bubble already present, to 'nucleate' other 
clusters/stripes/bubbles throughout the remainder of the 
system, (ii) Set the initial input density profile to be 



p(r) =p + x (r), 



(19) 



where x(r) is a random noise field with the property 
| x(r) | <C p. This is the "quick and dirty" way of estab- 
lishing the fluid phase behaviour for a particular p. (iii) 
Take the density profile obtained previously for a neigh- 
bouring state point and use it as an initial guess for the 
density profile at the required state point. Approaches 
(i) and (ii) proved to be good for an initial assessment of 
the phase behaviour. However, none of these approaches 
guarantees that one finds the density profile correspond- 
ing to the global minimum of the grand potential. To 
find the true global minimum one must follow approach 
(iii) above for all the possible candidate structures and 
then take the resulting density profiles (which each corre- 
spond to local minima in the free energy landscape) and 
substitute the density profiles into the grand potential 
functional ^ and evaluate the grand potential f2 corre- 
sponding to each phase. The phase with the lowest grand 
potential is the true equilibrium phase. Note also that 
since we were solving for the density profile on a finite 
grid of size L x x L y , with periodic boundary conditions, 
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FIG. 2: (Color online) Typical density profiles in each of 
the microphase separated phases. These profiles are for the 
temperature ksT/e = 0.15 and the average densities are 
pa 2 = 0.2, 0.35 and 0.5 in the cluster, stripe and bubble 
phases respectively. The x- and y-axes are in units of a. 



one must also minimise the grand potential with respect 
to L x and L y as well as with respect to variations in 
p(r), since otherwise one is forcing the period of the fluid 
modulations to be L x /n in the x-direction and L y /m in 
the y-direction, where n and m are integers [27 . This 
process is not as arduous as it might first appear. This is 
due to the fact that, as noted in the previous section, the 
length scale determining the size of the periodic struc- 
tures is 2ir/k c , which for the present DFT does not vary 
as a function of density or temperature. This means that 
having found the values of L x and L y which minimise f2 
for a certain state point, the values of L x and L y which 
minimise f2 for a different point in the phase diagram are 



very close to those at the first state point. 

In order to locate the phase transitions between the 
various phases and calculate the densities of the coexist- 
ing phases, we performed scans of the chemical potential 
p for fixed temperature T. Along these scans we recorded 
the pressure P (obtained from the grand potential, since 
fl = —PA, where A is the system area) and the aver- 
age density of the system p. We repeat this for all the 
phases displayed by the system. Since phase coexistence 
occurs between points with equal (T, P, pi) one can then 
read off the coexisting state points by plotting P versus 
pt, for these isothermal scans and noting the intersection 
points between the different curves. The resulting phase 
diagram is displayed in Fig. [T] In addition to the homo- 
geneous vapour and liquid phases we find that the DFT 
predicts that the system displays three inhomogeneous 
fluid phases. Examples of density profiles obtained for 
each of these three phases are displayed in Fig. [2j At 
lower densities, the system forms a regular array of clus- 
ters. The density within the clusters is close to that of the 
uniform liquid for the same temperature and the density 
in the voids is low - close to that of the uniform vapour 
for the same temperature. At intermediate densities, the 
system forms an array of parallel stripes. Again, the par- 
ticle density within the stripes is that of a liquid and the 
local density between the stripes is that of the vapour. 
At higher densities still, the system forms a regular array 
of bubbles, i.e. within the bubbles the density of parti- 
cles is low, but in between the bubbles the local density 
is that of the liquid. 

All of the phase transitions are predicted by the present 
DFT to be first order, except for the transition from the 
uniform fluid to the stripe phase. This phase transition is 
predicted to be second order (continuous) . If one starts 
in the uniform fluid phase at a state point directly above 
the stripe phase in temperature, and then decreases the 
temperature, one finds that the system is uniform right 
up to the A-line, which is itself the phase transition line. 
The relevant order parameter for the phase transition is 
the amplitude of the stripe modulations A [H] . Outside 
the A-line, and on the A-line itself, the amplitude .4=0. 
On decreasing the temperature below T\, the tempera- 
ture on the A-line, one finds that the amplitude of the 
modulations in the stripe phase grows continuously with 
decreasing temperature |21j . We will discuss this issue in 
further detail in Sec. IVII 



V. CONFINED FLUID PROPERTIES 

In this section we apply the DFT Q to investigate 
the behaviour of the fluid under confinement. In par- 
ticular, we examine the case when the fluid is con- 
fined between two parallel hard walls, where the exter- 
nal potential varies in only one Cartesian direction and 
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V ext (r) = V ext (x), where 




x < 
< x < L 
x > L 



(20) 



This situation was studied in Ref. [13] by Imperio and 
Reatto using Monte Carlo computer simulations. As 
in the previous section, we calculate the equilibrium 
fluid density profile p(r) by discretising it on a two- 
dimensional Cartesian grid and then solving Eq. ( 18 1 via 



a simple iterative numerical procedure, starting from an 
initial guess for the density profile. This is done using pe- 
riodic boundary conditions in the Cartesian y-direction 
(the direction parallel to the walls). 

For a given temperature T, we investigate the various 
morphologies that the fluid density profile displays as the 
distance between the two walls, L, is increased. These 
calculations could be performed for fixed chemical po- 
tential, /i. However, the Monte-Carlo simulation results 
in Ref. [T3] were performed in the NVT ensemble (fixed 
number of particles N, volume V and temperature T), 
keeping the average particle density between the walls 
fixed, as L was varied. In order to compare with these 
results, we perform our calculations at fixed average den- 
sity p. Thus, for each value of L, the chemical potential 
is chosen so as to achieve this target surface density. This 
is done by renormalising the density profile to the desired 
value at each step in the iterative routine for calculating 
the density profile. This method means that one does 
not have to know a priori the precise value of p required 
to achieve the desired average density. 

For a given L, one finds that there can be several dif- 
ferent density profiles which correspond to minima of the 
free energy - such as, for example, one density profile ex- 
hibiting stripes parallel to the walls and another with 
stripes lying perpendicular to the walls. As in the previ- 
ous section, in order to find the various candidate density 
profiles we implemented the following three different ap- 
proaches: (i) choose an initial guess for the input density 
profile with a cluster/stripe/bubble already present, to 
nucleate other clusters/stripes/bubbles throughout the 
remainder of the system, (ii) Starting from the noisy 
initial density profile in Eq. (19 1, (iii) Take the den- 



sity profile obtained previously for a neighbouring state 
point or value of L as the initial guess. The approach 
taken was to use methods (i) and (ii) to generate the set 
of candidate structures and then we used method (iii), 
to calculate the Helmholtz free energy via Eq. Q, as a 
function of L. 

In Fig. [3] we display the Helmholtz free energy per unit 
of available area between the walls A, as a function of 
the distance between the walls L (upper figure), together 
with a sequence of representative density profiles (lower 
figure). The temperature is ksT/e = 0.2 and the average 
density is pa 2 = 0.2, which as can be seen from Fig. [l] 
corresponds in the bulk to the cluster phase. For small 
values of L < 2ir/k c , the fluid density profile does not 




FIG. 3: (Color online) Top: the Helmholtz free energy per 
unit area, T/A, as a function of slit width, L, for the case 
when the average fluid density in the slit is pa 2 = 0.2 and 
the temperature ksT/e = 0.2. Bottom: a series of density 
profiles for increasing L, labelled (a) - (e) which correspond 
to the different free energy curves in the upper figure, which 
are correspondingly labelled. 



vary in the y-direction (parallel to the walls) . The branch 
of the free energy corresponding to this case is labelled (a) 
in the upper part of Fig. [3] and a typical density profile for 
this branch of the free energy is correspondingly labelled 
in the lower part of Fig. [3] As L is increased, we find a 
new minimum of the free energy, which corresponds to a 
density profile exhibiting a single line of clusters - see the 
density profile labelled (b) in the lower part of Fig. [3] and 
the correspondingly labelled free energy curve above. On 
further increasing L, one finds a sequence different free 
energy branches, with the density profile along each or 
these, having an additional line of clusters between the 
walls. Free energy curve (c) corresponds to two lines 
of clusters lying parallel to the walls, (d) to three lines 
of clusters, (e) to four lines of clusters and (f) to five 
lines of clusters (the corresponding density profile is not 
displayed in this case). 

In Fig. [4] we display the Helmholtz free energy per unit 
area and a number of typical density profiles for the case 
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FIG. 4: (Color online) Same as Fig. [3] except here the average 
density in the slit is per 2 = 0.3. 
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FIG. 5: (Color online) Some of the metastable density pro- 
files that were found for the case when fcsT/e = 0.2 and the 
average density pa 2 = 0.3. 



when ksT/e = 0.2 and pa 2 — 0.3. In bulk, this state 
point lies within the stripe phase and so for most values 
of L the equilibrium density profile consists of stripes that 
lie parallel to the walls. However, for values of L where 
the free energies for a density profile with n and for a 
density profile with (n + 1) stripes arc equal (i.e. when L 
is incommensurate with the stripes), then one observes 
that the density profiles exhibiting clusters have a lower 
free energy and are the equilibrium configurations. Thus, 
even though for this average density the bulk phase is the 
stripe phase, due to the fact that this state point is not 
too far in the phase diagram from the cluster phase, one 
sees the influence of this proximity in the confined fluid 
density profiles. For small values of L, the fluid density 
profiles do not vary in the direction parallel to the walls 
- the branch of the free energy corresponding to this is 
labelled (a) in the upper part of Fig. [3] On increasing 
L we find (b), corresponding to short stripes lying per- 
pendicular to the wall, followed by (c) where we find two 
stripes lying parallel to the walls. Increasing L, we find 
(d) corresponding to two parallel stripes with a line of 
clusters in between. Further increasing L we find (e), 
corresponding to three stripes lying parallel to the walls. 
As L is increased even further we see an increasingly com- 
plex series of transitions between configurations with dif- 
ferent morphologies - see the magnification in the inset 
of the upper panel of Fig. |4j Note that the density pro- 
files displayed in Fig. [4] are the ones that correspond to 
equilibrium density profiles (i.e. to global minima of the 
free energy, for some value of L). However, we also find 
a large number of other density profiles with higher free 
energies; see for example some of these displayed in Fig. 
[5] These metastable configurations correspond to local 
minima in the free energy that are not global minima. 

In Fig. [6] we display the free energy per unit area to- 
gether with a number of representative density profiles for 
the case when ksT/e = 0.2 and pa 2 = 0.4. In bulk, this 
corresponds to a state point in the bubble phase, lying 
very close to the transition to the stripe phase. How- 
ever, due to the influence of the confining walls, we find 
that for most values of L investigated here, the equilib- 
rium density profiles consist of stripes lying parallel to 
the walls. Thus, the confining walls exert a strong in- 
fluence on the symmetry of the fluid density profiles, for 
this state point. Note that the 'ladder' density profile 
labelled (x) corresponds to a configuration that does not 
correspond to the global minimum to the free energy It 
is included in Fig. [6] because for some values of L it has a 
free energy that is very close to that of the true minima 
profiles, (e) and (g) - see the magnification in the inset 
of the upper panel of Fig. [6j 

In Fig. [7] we display results for the case when fc^T/e = 
0.2 and pa 2 = 0.5. In bulk, this corresponds to a state 
point in the bubble phase, lying well removed from the 
transition to the stripe phase. However, due to the influ- 
ence of the confining walls, for small values of L we find 
that the equilibrium density profiles consist of stripes ly- 
ing parallel to the walls. Note also that the transition 
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FIG. 6: (Color online) Same as Fig. [3] except here the average 
density in the slit is pa 2 = 0.4. 



from (a), having one stripe, to (b), having two stripes, 
appears to be a continuous transition. The metastable 
free energy branch labelled (x) in Fig. [7] is that corre- 
sponding to three, four and then five (as L is increased) 
stripes lying parallel to the walls. The transition from 
three to four stripes and also from four to five stripes 
appear to be continuous transitions. Interestingly, this 
is not the case for the transition between two and three 
stripes. For L/a > 21 the equilibrium density profiles 
contain lines of bubbles lying parallel to the walls. 

In Ref. [T3] Imperio and Reatto used Monte-Carlo 
computer simulations to investigate the behaviour of the 
present model fluid under confinement. Due to the fact 
that there is only qualitative agreement between the bulk 
phase diagram predicted by the present DFT and that 
obtained from simulations, it is hard to compare results 
directly. In Fig. 1 of their paper, Imperio and Reatto 
display the potential energy as a function of L and also 
display a sequence of simulation snapshots for the case 
when pa 2 — 0.4 and the temperature T ~ 0.8T max , where 
T max ~ 0.13e/fcs is the maximum temperature at which 
modulated phases are observed to occur in the simula- 
tion results (for the present DFT, T max ~ 0.26e/fcs; see 
Fig. [I]). The sequence of structures displayed in Fig. 1 
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FIG. 7: (Color online) Same as Fig. [3] except here the average 
density in the slit is pa 2 = 0.5. 



of Ref. [13] do not correspond exactly with any of the 
results displayed here in Figs. [3] [4] [6] or [7j However, for 
values of L < 30a, the sequence of phases observed in 
the simulation results displayed in Fig. 1 of Ref. [13] are 
the same as those displayed in Fig. [4] From this com- 
parison, and also from comparing with the results in Fig. 
[6] which is a state point with the same average density 
as the results in Fig. 1 of Ref. [13], we conclude that the 
results in Ref. [13 effectively correspond to a state point 
somewhere in-between the two state points for which re- 
sults are displayed in Figs. [4] and [6j Since a detailed 
comparison between our results and the simulation re- 
sults of Ref. [13] is not possible, we are only able to draw 
the general conclusion that the present theory appears to 
be qualitatively reliable for describing the confined fluid 
properties of the present system and seems to at least be 
able to describe some of the sequences of structures that 
are observed in the confined fluid. 



VI. DISCUSSION AND CONCLUSIONS 

In this paper we have presented a simple DFT for 
describing the bulk and inhomogeneous fluid phase be- 
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haviour of a two-dimensional fluid with competing inter- 
actions. In Refs. [lOj [IT] the authors used Monte Carlo 
computer simulations in the NVT ensemble to study the 
structure and phase behaviour displayed by the present 
model fluid. The DFT results presented here are gener- 
ally in good qualitative agreement with these simulation 
results, as regards the topology of the phase diagram and 
the structure of the fluid in the various uniform and in- 
homogeneous phases. However, the present DFT is not 
able to describe some properties of the fluid. We discuss 
these weaknesses now. 

The present theory over estimates by a factor of 
two the maximum temperature at which microphase- 
separated phases are expected to occur. This is largely 
due to the fact that in the present theory we retain for 
all values of r the form of the potential w(r), given in 
Eq. Q, in the mean field contribution to the free energy 
in Eq. Q. The value of w(r) for r < a (within the hard 
core overlap distance) should not contribute to the free 
energy. That the value of w(r) for r < a does give a 
contribution to the free energy, is a consequence of the 
mean field approximation made in constructing the free 
energy ^ - i.e. this equation predicts that the free en- 
ergy of the uniform fluid contains a term proportional to 
a = — J drw(r), the integral over w(r) for all values of r. 
A better approximation would be to truncate the poten- 
tial w(r) for r < a. This would give much better agree- 
ment with the simulation phase diagram 20, 21 , modify- 
ing the value of a = —w(0) and effectively just re-scaling 
the temperature axis, without qualitatively changing any 
of the results of this paper. We chose not to do this, how- 
ever, due to that fact that the analysis given in Sec. Ill 
for quantities such as S(k) and k c , is much less straight- 
forward if one truncates w(r), because one then finds that 
w(k) no longer has the simple form given in Eq. ( 11 ). 



More substantial differences between the present the- 
oretical predictions and the simulation results are as fol- 
lows: In Ref. jTTj the authors show that for temperatures 
T < T max , on increasing the density there is a transition 
from the vapour phase to a cluster fluid phase - i.e. to a 
disordered system with no long range ordering of the clus- 
ters. On further increasing the density they then found 
a further transition to a cluster phase with crystalline 
ordering of the clusters, in line with the predictions from 
the model presented here. The present theory predicts 
only a single cluster phase exhibiting crystalline ordering. 
We believe this is due to the fact our DFT is a simple 
mean field theory, which neglects certain fluctuation con- 
tributions to the free energy [25]. These neglected con- 
tributions must play an important role in the transition 
from the disordered cluster fluid phase to the ordered 
cluster phase. 

A second difference between results from our DFT and 
the simulation results in Refs. [TO] [TT] are that for tem- 
peratures T < 0.6T max , the system starts to display crys- 
talline ordering of the particles within the clusters and 
stripes. Owing to the fact that we have made a local den- 
sity approximation in the free energy functional ^ for 
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FIG. 8: (Color online) Density profiles calculated for 
fcsT/e = 0.2. For the upper profile, the average density is 
pa 2 — 0.3, a state point in the stripe phase. The lower profile 
is for pa 2 = 0.45, a point in the bubble phase. These profiles 
were obtained by quenching from a homogeneous system with 
a small amount of random noise added to the density profiles 
to nucleate the structuring - see Eq. ( 19 1. In the stripe phase 
this results in a density profile containing many defects. Sim- 
ilarly, in the lower density profile this results a number of 
domains forming in this (periodic) system. The true mini- 
mum in the grand potential corresponds to periodic arrays 
with no defects. 



the hard-disk contribution to the free energy, our theory 
is unable to describe this effect. To describe such freezing 
effects one must implement a much more sophisticated 
(non-local) hard-disk reference free energy functional. 

We return now to the issue raised at the end Sec.[lV] re- 
garding the order of the phase transitions. In the present 
DFT, all the phase transitions between the equilibrium 
phases are first order with the exception of the transition 
from the uniform fluid to the stripe phase, which is pre- 
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dieted to be a continuous phase transition. To further 
understand the nature of the phase transition between 
the uniform fluid and the stripe phase, we may follow 
the approach used in Ref. [2T]. Starting from the DFT 
p} . and assuming that the density profile is of the form 



p(z) = p + Asin(kz), 



(21) 



the free energy may be expanded in powers of the order 
parameter A in the stripe phase. One finds that the 
Helmholtz free energy then takes the form 

F[p(z)} = T[p\ + b 2 (p, k)A 2 + b 4 (p)A 4 + 0(A 6 ). (22) 



The coefficient b 2 oc a(r/) + p(3w{k) [c.f. Eq. ( 10 )] and 
the coefficients b n > 0, where n > 4. The value of the 



amplitude A which minimises the free energy ( 22 1 is the 
solution to the equation dT/dA = 0. Outside the A- 
line, the coefficient 62 > for all values of k, so that 
the minimum of the free energy is when A = 0, i.e. for 
the uniform fluid. On the A-line itself, b2(p\,k c ) = 
and within the lambda line 62 {p, k c ) < 0, indicating that 
within the A-linc the modulated fluid has a free energy 
that is lower than that of the uniform fluid. At higher 
temperatures this transition from the uniform fluid to the 
stripe phase is a second order phase transition, and the 
transition line is the A-line itself. Both the amplitude A 
and the density p vary continuously across the transition 
line (although the first derivative of A(p, T), does not 
|21|). On descending to lower temperatures, one finds 
that the transition between these two phases ceases to 
be continuous and becomes first order, due to the density 
and k dependence of the coefficient b 2 (p, k), and there are 
two tricritical points connected to one-another by the A- 
line [3T]. For temperatures below these tricritical points, 
both the density p and the amplitude A vary discontinu- 
ously at the phase transition. These tricritical points are 
not displayed in Fig. [TJ because these points are located 
in the regions of the phase diagram where the cluster and 
bubble phases are the equilibrium phases. Thus, in the 
present two-dimentional system, these tricritical points 
are not accessible, in contrast, it seems, to the case in 
three-dimensions [9| I21j. 

In the simulation studies of Imperio and Reatto [10, 
ITT] , the phase transitions from the uniform fluid to the 
modulated fluid phases are located by calculating the ex- 
cess specific heat C% x as a function of temperature. This 
quantity displays a peak and the phase transition is iden- 
tified with this peak. They find that for low and high 
densities p, the peak is well pronounced. However, for 
intermediate densities around pa 2 ~ 0.3, they find that 
the peak is greatly reduced in amplitude, in comparison 
with the height of the peak at higher or lower densities 
|llj . Comparing these findings with the present DFT 
results, we see that where the DFT predicts the phase 
transition from the uniform fluid to the stripe phase to 
be second order, coincides with where there is almost no 
peak in C% x , and densities where the DFT predicts the 
transition from the uniform fluid to the modulated fluid 



to be first order, the simulation results show that there 
is a pronounced peak in C% x . In their analysis of the 
system size scaling of C% x , Imperio and Reatto conclude 
that their results for the densities pa 2 — 0.05, 0.15 and 
0.4 point to the possibility of the phase transition from 
the modulated fluid phase to the uniform fluid phase be- 
ing a Kosterlitz-Thouless transition |llj . In order to con- 
firm this theoretically, one must go beyond the present 
mean field DFT theory. One possible route to do this is 
perhaps by following the approach of Ref. [25] . 

In the present two dimensional system, fluctuations 
are not only important in determining the nature of the 
phase transitions, but are also important in determining 
the structure of the bulk modulated phases themselves. 
Let us consider in particular the stripe phase, although 
the discussion that follows may also have implications 
for the cluster and bubble phases. The present DFT 
predicts that the equilibrium configuration in the stripe 
phase consists of an array of perfectly parallel stripes, 
displaying long range order. However, as we have noted 
already, the present DFT is a mean field theory and ne- 
glects certain fluctuation contributions to the free energy. 
Thus, although the minimum of the free energy corre- 
sponds to perfectly parallel stripes, in reality one should 
find that fluctuations about this minimum will destroy 
the long range order - see for example the discussion in 
Refs. [55J[3U]. In fact, when one seeks for the minimum 
of the free energy, using as initial guess in the numeri- 
cal procedure the noisy density profile given in ( 19 1, one 
finds that as the system size L x x L y increases, there 
are an increasingly large number of metastable minima 
in the free energy. The density profiles for these minima 
may contain many defects - see for example the results 
displayed in Fig. [8j These density profiles are calculated 
for state points in the stripe phase (upper figure) and in 
the bubble phase (lower figure). Whilst the DFT shows 
that such configurations corresponds to a minimum in the 
free energy, the theory does not reveal the height of the 
free energy barrier between this configuration and any 
other neighbouring configurations. In order to address 
the influence of fluctuations and to tackle the question 
of the size of the barriers between free energy minima, 
we propose to extend the present theory in the following 
manner: If one assumes over-damped stochastic equa- 
tions of motion for the colloidal particles, then following 
the approach of Ref. |31j one may develop a dynamical 
DFT, which predicts the following stochastic equation 
of motion for the time dependence of the coarse grained 
fluid density profile: 



dp(r,t) 
dt 



= V 



p(r,t)V 



+ V-y/p(r,t)£(r,t) 



5p(r,t) 

_ (23) 

where T[p] is the free energy functional in Eq. Q and 
£(r, t) is Gaussian random noise field. Such an approach 
may allow one to determine the influence of fluctuations 
in the present system and perhaps also to address the 
question of whether the fluid effectively becomes trapped 
in configurations such as those displayed in Fig. [8] as the 
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temperature is decreased - i.e. to determine if and where 
the glass transition line is in the present system [32j [33] . 
We plan to pursue this line of investigation in the future. 

The sequence of structures displayed by the present 
model fluid (clusters to stripes to bubbles), has also 
been observed in the configurations of a subcritical con- 
strained lattice gas (Ising model) [M] and a fluid of rod- 
like particles |35] - both of these are two-dimensional 
systems. Similar behaviour has also been observed in 
the three-dimensional Lennard- Jones fluid |36j . None 
of these models have long-ranged repulsive interactions 
between the particles, in contrast to the present sys- 
tem. This sequence of structures is only observed when 
these systems are constrained within a finite-sized box 
with periodic boundary conditions and with the num- 
ber of particles fixed, so that the average density in the 
box is between that of the coexisting vapour and liquid 
phases (34j [35j [36j ■ The cluster/stripe/bubble structures 
in these systems arise due to the constraint that the num- 
ber of particles is fixed and the structures are unstable 
with respect to fluctuations in the density. In contrast 
to this, in the present system the configurations that we 
observe are stable equilibrium structures and do survive 
in an open (grand-cannonical) system. The long range 
repulsion between the particles stabilises the structures, 
in just the same way as the long-range repulsion stabilises 



the structures exhibited by three-dimensional fluids with 
competing interactions [9]. A further difference is that 
in systems without the long range repulsive interaction 
between the particles [331 I3HJ I3S], the size of the clus- 
ters/stripes/bubbles depend on the size of the box they 
are in. However, in the present system, it is the com- 
petition between the short range attraction and the long 
range repulsion that determines the size of the structures 
and as long as they are within a sufficiently large box, the 
size of the clusters/stripes/bubbles is independent of the 
system size. 

To conclude, we remind the reader that there are 
many important applications for systems having col- 
loidal nanoparticles that are confined to a fluid inter- 
face, ranging from optical devices to stabilising emulsions 
[2\ . Whilst the present and other studies have gone some 
way towards understanding the self-assembly of such two- 
dimensional fluids, there arc still several open questions 
to be addressed [2J. 
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